Taylor dispersion with absorbing boundaries: A Stochastic Approach 
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We describe how to solve the problem of Taylor dispersion in the presence of absorbing bound- 
aries using an exact stochastic formulation. In addition to providing a clear stochastic picture of 
Taylor dispersion, our method leads to closed-form expressions for all the moments of the convective 
displacement of the dispersing particles in terms of the transverse diffusion eigenmodes. We also 
find that the cumulants grow asymptotically linearly with time, ensuring a Gaussian distribution 
in the long-time limit. As a demonstration of the technique, the first two longitudinal cumulants 
(yielding respectively the effective velocity and the Taylor diffusion constant) as well as the skew- 
ness (a measure of the deviation from normality) are calculated for fluid flow in the parallel plate 
geometry. We find that the effective velocity and the skewness (which is negative in this case) are 
enhanced while Taylor dispersion is suppressed due to absorption at the boundary. 
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In a laminarly flowing fluid with flow velocity varying 
perpendicular to the flow direction, suspended particles 
will jump randomly across fast and slow streamlines due 
to transverse diffusion. This fluctuating convectional ve- 
locity arising from a random sampling of various stream- 
lines disperses these particles along the direction of flow 
in a well-studied process known as Taylor dispersion[l- 
4]. The particles drift along the flow with an average 
velocity v e and disperse linearly with time t along the 
flow, i. e. the r.m.s displacement in that direction is 
\J (-Dtayior + D)t where the convection-induced Taylor dis- 
persion coefficient D taylor is found to scale as £ 2 v^/D, D 
being the molecular diffusion coefficient and I being the 
transverse dimension. Taylor dispersion has turned out 
to be important both fundamentally in hydrodynamics 
and in its applications to diverse fields ranging from bi- 
ological perfusion, chemical reactors and lab-on-chips to 
soil remediation and oil recovery [4- II]. 

In spite of having been around for many decades, Tay- 
lor dispersion in the presence of absorption, even in sim- 
ple geometries (such as Poiseuille flow in tubes), has 
proven to be a mathematically challenging problem [ 12- 
16]. As Taylor dispersion under Poiseuille flow in tubes 
is now routinely used to measure molecular diffusion 
coefficients [17], it is important to compute the errors due 
to absorption at the walls. The same holds true for other 
applications outlined in the previous paragraph. The 
methods used originally [1-3] worked only in the absence 
of absorption. This was pointed out in the most 'modern' 
treatment available — by Gill, Lungu and Moffatt[15, 16] 
dating back to over 40 years ago. They showed that the 
effective velocity is enhanced while Taylor dispersion is 
suppressed due to absorption at the boundary. However, 
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even in their scheme, obtaining higher than second mo- 
ments and thus estimating deviation from Gaussianity of 
the particle concentration profile along the direction of 
flow, becomes prohibitively tedious and hence was not 
addressed. 

Our approach goes a step further than this cornerstone 
work. We provide the correct stochastic picture of the 
underlying process. When decay is present, the usual 
probabilistic techniques [12, 13] fail. We formulate, for 
the first time, the correct technique — ultimately pro- 
viding a visually appealing and computationally simple 
stochastic method to find all moments / cumulants of the 
longitudinal (i. e. along the convective flow direction) dis- 
placement for uniform (linear) laminar flow. We find that 
the cumulants grow asymptotically linearly with time en- 
suring a Gaussian distribution in the long-time limit, as is 
characteristic of drifting random walkers on a lattice[18]. 

Consider diffusing particles being carried along by a 
fluid flowing unidirectionally in the x-direction through 
a uniform cross-section A. The particle position in the 
transverse section A will be denoted by y, while the com- 
plete position will be specified by x = (x,y). The con- 
centration of the particles, N(x,t), evolves according to 
the convective diffusion equation inside the fluid volume 
V: 

dN<K Q t '^ = D VlN(x, t) - v(x).VN(x, t) xeV (1) 

where D is the molecular diffusion constant of the par- 
ticles in the still fluid and v(x) = v(y)x is the x- 
independent convective velocity field of the fluid. The 
presence of absorbing walls is taken into account by the 
boundary condition 

De.V x N(x,t)+pN(x,t) = 0, x G dV (2) 

where p is a measure of the surface absorption rate that 
varies from zero (reflecting walls) to infinity (perfect ab- 
sorption), while e(x) is the normal to the surface dV at 
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x. Since the cross-section is uniform, e(x).x = 0. This 
enables us to replace e.V x by e.Vj, in (2). 

By integrating over x in (1) and (2), a diffusion equa- 
tion for n(y,t) — J dxN(x,t) - the concentration of 
particles in the transverse section A - is obtained with 
absorbing boundary conditions (when N,d x N — ► as 
x — > ±oo): 

^=5V>(^) yei (3a) 

D e.V y n{y,t) = -pn{y,t) y e dA (3b) 

This yields the y-statistics of the dispersing particles that 
we will use in the formalism below involving the equa- 
tions (4) and (5). 

We now proceed to find the moments of x(t). We shall 
not show explicitly the noise that is responsible for the 
usual molecular diffusion along the x direction and con- 
centrate on the dominant convective noise [21] encoded in 
the stochastic equation that gives rise to the convective 
part in (1). This equation is 

^ = v(y(t)) => x{t) = J* dt'v{y{t')) (4) 

Intrinsic diffusive noise along the x-direction may be in- 
corporated easily — in particular, the mean will be un- 
affected while -Dtayior will be additively augmented by D. 
We are considering the case ir(0) = here. To proceed 
further, we introduce the Green's function/propagator 
G(y, t\i?) - the t > solution to (3) when n(y, 0) = 
5{y — y 1 ). It's most important use is encoded in the 



formula[19] 

n(y, t) = J dy f G(y, t - i'|y>(^, *'). (5) 

for t > t' . In other words, it gives us the fraction of par- 
ticles starting out from an arbitrary initial state (y',t') 
and surviving till a specific final state (y, t) . 

To illustrate the technique of obtaining the moments 
of x(t), we start with the first one - the mean. It is ob- 
tained by averaging (4) over the particles that survive 
the absorbing walls till time t. Following (5), the parti- 
cle density at (yV) is given by / G(y,t'\y i )n(y i ,0)dy l . 
However, in the remaining time (t — t') only a fraction 
G(yf,t — t'\y) of these will survive till a final state (yf,t). 
Thus, the number density of particles that pass through 
(y 1 ',£') and also survive till t is 

v{y,t'\t) - f dy f f d yi G{y f ,t - t'lyjG&t'lyAniy^O) 

J A J A 

Multiplying this density by the corresponding displace- 
ment v(y)dt', averaging over all possible intermediate po- 
sitions and finally adding up all these averaged displace- 
ments, we find that 

<*(*)> = W{t)f dt 'Sfi "GTWI*) (7) 

where the (normalizing) denominator N(t) is equal to the 
total number of particles that survive till t. Similarly, 
squaring (4) and averaging over the surviving particles 
one can show that the second moment is given by 



{{x(t)f) 



v(yi,ti\y2,h\t) 



N(t) 



A® A 



/ dhdt 2 / dy 1 dy2v(y 1 )v(y 2 )i>(yi,t 1 \y2,t2\t) 

J0<t 1 <t 2 <t J A® A 

dyidy f G(yf,t - t 2 \y 2 )G(y 2 ,t 2 - ti\yi)G(yi,ti\yi)n(yi,0) 



(8a) 
(8b) 



The density of walkers who survive till t after passing 
through (j7i,ii) and (y 2 ,t 2 ) is given by v{yi,h\y2,t 2 \t), 
analogous to (6). All higher moments may be obtained 
by repeating the arguments leading to the above results. 



The preceeding arguments thus provide the correct ex- 
position of the stochastic technique to be used when ab- 
sorption is also present. This will be the mainstay of the 
remaining results to be derived in this paper. We note 
that the formulation is exact in the framework we con- 
sider — viz. the statistics obeyed by y is independent of 
x. 



We can solve (3) by expanding in terms of diffusion 



eigenmodes[19]: 



i(y,t) =^n k ^ k {y)e Afct 



fc=0 



where no ^ (necessarily) and 



(DV 2 y + x k )Mv) = o 

(De.V y + p)My) = 
(A < Ai < A 2 . . .) 



dy ipj{y)i>k{y) = 5 jk 



ye A 
yedA 



(9) 



(10a) 
(10b) 

(10c) 
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The propagator is given by [19] 

oo 



(ii) 



71=0 



The lowest eigenvalue Ao for this problem becomes non- 
zero for p ^ — this is the main difference between the 
absorbing and non-absorbing cases. At long times all the 
higher modes decay with inverse mean lifetimes of 



A,, 



n 2 /r 



(12) 



where t is the characteristic relaxation time-scale of the 
transverse diffusion process and scales as £ 2 /D where £ is 
the characteristic length-scale of the cross-section A. In 
the long-time limit Orwe can thus work with only the 
lowest mode in (9) and the errors will go as 0(e~*/ T ). 

We consider this long time limit in what follows. We 
find that the moments can be conveniently evaluated in 
terms of simple 'matrix elements' 



\ dy ipj(y)v(y)ip k (y) (13) 

J A 



Using (11) and (10) in (7) and keeping only the dominant 
terms as outlined in the previous paragraph we obtain the 
following expression for the first moment in the long time 
limit: 



(x(t)) ^^<^t + 0(e- 



) 



(14) 



This result is very counter-intuitive since one naively ex- 
pects, at late enough times, the effective velocity v e to 
be the same as the asymptotic instantaneous mean ve- 
locity that is proportional to j A v n cx J A vip — since 
n(y, t — > oo) ~ ipo(y)- However, absorption at the walls 
modifies the probability that a particle starting out at 
y survives after a long time (i. e. 3> r) from being y- 
independent to being proportional to ipo(y) — this pro- 
vides the additional power of Vo in (14) (recall that 
vqo =fvi/>$). 

Analogously, this time using (8) instead of (7) above 
we calculate the variance in x to be 



<(<fc(*)r> 



2\ *>3r 



t X 2 



fc=i 



u k0 



+0(e- t/T ) (15) 



Otaylor 



while the third cumulant simplifies to: 

OO OO / r- \ 

\3\ *>4t , \ - \ - VOjjVjk - OjkVQOjVkO 



(WW) 



~ik~i _ x o)(^k - A ) 



(16) 



showing that the skewness, which is a measure of the 
deviation from Gaussianity, dies out with time: 



((faffl) 3 ) 
((,5x(t)) 2 ) 3 /2 



(17) 



Using the same arguments, higher cumulants may be sim- 
plified to forms analogous to equations (14) or (15) in the 
long time limit. Because of the straightforward algebra 
involved, computer algebra systems (we used Mathemat- 
ica) may be programmed to do this correctly and elimi- 
nate corrections 0(e~*/ T ). 

We have analytically computed uptil the fifth cumulant 
this way and found that they are all proportional to time 
in the long-time limit. We assert that this holds to all 
orders. This implies that the p.d.f of the longitudinal 
displacement tends to a Gaussian with deviations that 
die as 1/t — similar to the case of reflecting boundary 
conditions [3]. A similar behaviour is also exhibited by 
drifting random walkers on a lattice[18]. 




v x {y) =v (l-y 2 /e 2 ) 



y = +e 



y = fluid 



V = 



FIG. 1: Taylor dispersion between parallel plates 

Finally let us compute a specific example - Taylor dis- 
persion of particles suspended in fluid flowing unidircc- 
tionally between two infinite parallel plates (see fig. 1). 
Absorption at the walls is characterised by the dimen- 
sionless parameter a — pl/D (p is defined in (2)). For 
convenience, we have used the dimensionless coordinate 
1] = y/£ below. The eigenf unctions (10) for this case turn 
out to be: 



(18a) 
(18b) 




n =0 



(19a) 



™ + ?l-JL(?EL) + o(—X 

2 n7T nil \niT J \n-K J 



n =1 



For a — > oo, 



(n + 1)tt 



0,1,... 



(19b) 



(20) 



Finally, the previously introduced eigenvalues A„ (sec 
(10)) are obtained from the relation: 



X n =D 



(21) 
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Using our formulae (14) and (15) we have calculated 
the effective velocity and the Taylor dispersion constant 
to be: 



Ve = ir (' • j^" • c,: ""' : ) ; 

2v a ( 3 



(22a) 



-D tiivlor 



' : = 945~D I 1 " Y5 a + ° {a ))' a<<1 



= 14% of value at a = 0, when a — > oo (22b) 

When there is no absorption (a = 0), w e is equal to the 
asymptotic mean instantaneous velocity of the particles 
since ip is constant. However, since the slower particles 
near the walls are also the ones that are preferentially 
absorbed, v e increases with a as only surviving parti- 
cles contribute to the moments. Similarly, since more 
particles at the walls get removed for increasing a, the 
dispersion in the convection field as seen by the particles 
goes down and so does D tayloI which arises due to it. Our 
results (22) agree with those derived in [16] by a more 
tedious method. One may reproduce (15) from [16] with 
some effort. However, as already noted before, it is very 
difficult to obtain simple formulas for the higher moments 
using the method presented there. 

Finally, the skewness, which is a measure of the devi- 
ation from gaussianity, is found to be 



showing that it is negative and enhanced due to absorp- 
tion at the walls. The zero-absorption case of this result 
(23) has been independently verified using the method 
presented in [3]. 

In conclusion, we have a picturesque and computation- 
ally simple solution to a problem of general interest. For 
example, we may now estimate the error in measuring 
D using Taylor dispersion [17]. The value of Ao and con- 
sequently a may be obtained from the observed expo- 
nential decay of the particle numbers and subsequently 
used in the corrected formula (15) to estimate D from 
the experimentally-observed variance. Also note that the 
moments in the case of an arbitrary cross-section may be 
computed numerically to the necessary precision by in- 
volving a finite number of cigenmodes. This is possible 
because terms involving higher eigenmodes get progres- 
sively less important as may be easily seen from the forms 
of (14) and (15). 
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